Direct determination of the size of basins of attraction of jammed solids 
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We propose a free-energy based Monte-Carlo method to measure the volume of potential-energy 
basins in configuration space. Using this approach we can estimate the number of distinct potential- 
energy minima, even when this number is much too large to be sampled directly. We validate our 
approach by comparing our results with the direct enumeration of distinct jammed states in small 
packings of frictionless spheres. We find that the entropy of distinct packings is extensive and that 
the entropy of distinct hard-sphere packings must have a maximum as a function of packing fraction. 

PACS numbers: 61.43.-j,61.43.Bn,61.43.Fs 



When many equal-sized spheres are poured into a con- 
tainer, the spheres are unlikely to end up arranged in a 
periodic lattice. This observation reflects the fact that 
Sq, the entropy of distinct disordered packings that are 
mechanically stable, is very large compared to the corre- 
sponding entropy of distinct ordered packings. The fact 
that So is so large has important consequences for the 
disordered packings such as granular materials 

There is a natural connection between hard-sphere 
packings and glasses [1], whose potential energy land- 
scapes have many minima (inherent structures [5[), cor- 
responding to mechanically-stable states. These minima 
have been argued to be relevant for our understanding 
of the glass transition [J 0] . The number of such min- 
ima has been calculated from replica theory 0, 0, [l] ■ In 
calculating this number numerically, however, a protocol 
must always be used to generate energy minima. Typ- 
ical protocols produce states with probabilities that are 
not known; for example, when the entropy of minima 
is calculated from finite-temperature simulations ^-dH, 
one must assume that the temperature is low enough so 
that the system crosses no barriers. Similarly, when the 
entropy is calculated from algorithms that involve com- 
pression or dilation of the system [isl - fTsj , it may depend- 
even under ideal conditions-on the algorithm used. As 
a result, it is difficult to count the number of distinct 
mechanically-stable states, with states weighted equally. 

In this letter we report a general computational 
method to measure the volume of a basin of attraction 
associated with an arbitrary potential energy minimum. 
This is the key to calculating the entropy of distinct min- 
ima for soft spheres because there is a protocol that gen- 
erates minima weighted by their basin volumes In 
this "basin" protocol, states in configurational space are 
selected at random and each one is quenched to its near- 
est energy minimum By using this protocol and cor- 
recting for the weighting by calculating the basin vol- 
ume, we can obtain the unweighted entropy of distinct 
mechanically-stable states (packings) for soft spheres. Fi- 



nally, the analogous entropy for hard-sphere packings can 
be obtained from the density of soft-sphere packings at 
zero pressure. We find that there must be a maximum in 
the entropy of distinct hard-sphere packings, at least for 
small systems, in agreement with earlier results obtained 
by direct enumeration [isj . 

To explain our approach, we first define the volume of 
a basin for a packing of soft spheres as 



(1) 



vb = / dRGiR,Ro), 



where G{R, Rq) = 1 if, upon energy minimization, any 
point R in configuration space ends up at Rq, the position 
of the local potential energy minimum, and otherwise. 
The integral is over the whole configuration space. We 
view the (hyper) volume associated with a given basin 
as a partition function and hence compute its value by 
a suitable free-energy calculation method. Here, we will 
use the standard "Einstein" method [TB] and compute the 
basin free energy by comparing it to the free energy of 
a system confined near the minimum Rq by a harmonic 
potential with spring constant k. For arbitrary k, the 
canonical partition function of the system is: 



Q(fc) 



dRG{R, i?o)exp {-j3ku^ /2) 



(2) 



where u = \R—Rq\ is the distance between R and i?o, and 
/3 = (ksT)'^ with ks the Boltzmann constant. G{R, Rq) 
in Eq. @ can be rewritten as exp(— /3[/), where U = 
when R is in the basin, and oo otherwise. Obviously, 
Vb = Q(0). Without loss of generality, we choose (3 = 1. 

The free energy of this system is F{k) = —h\Q{k) and 
^^^^ = {v? /2)k, where {■■■)k denotes a canonical ensem- 
ble average at the spring constant k. This average can 
be sampled in a standard Monte-Carlo (MC) simulation. 
The change in free energy upon switching on a spring 
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constant km is 



F{km) = F{0) 



(3) 



where k„i is chosen sufficiently large that the confin- 
ing potential has no influence. In that case F{km) 
is known analytically and Eq. ([3]) allows us to com- 
pute F{0) and from that the volume of the basin, as 
Vb{Ro) = exp(— F(0)). In practice, we choose a max- 
imum kjn such that most (in our case > 90%) of the 
associated Gaussian distribution is within basin Rq. One 
then corrects the Einstein crystal result for the confining 
effect of the basin: F{km) ~ — ^In (27r/fc„i)— In/, where 
d is the dimension of space, N is the number of particles 
in the system, and / is the fraction of the associated 
Gaussian distribution within basin Rq. 

Given the basin volume, we can calculate the entropy 
of distinct mechanically-stable minima. We include in 
our analysis only energy minima that are mechanically- 
stable (jammed). The fraction of the total configuration 
space, Vtot .occupied by basins of jammed states, /j, is 



computed [J] by quenching randomly selected points in 
configuration space to the nearest energy minimum and 
calculating the fraction that end up in jammed states. 
The volume of configuration space at packing fraction (jj 
occupied by jammed basins is Vc{4>) = fji<t>)Vtot- 

As pointed out by Speedy in a slightly different con- 



text |17| . the total configuration space can be uniquely 
decomposed into distinct basins and hence its volume is 
simply the sum of the volumes of the constituent basins. 
Thus, 



^ = Oc I — ^ Vb^i I = f^c 

i=l \ 1=1 / 



(vb)- (4) 



By sampling the basin volume to obtain (wh), we can 
therefore compute Qc, the total number of distinct 
jammed states. 

Note that "computing the average basin volume" 
sounds simpler than it is because the probability to sam- 
ple a given basin is proportional to the basin volume 
itself. We correct for this bias by dividing by the basin 
volume. However, if a substantial fraction of all distinct 
basins together occupy a negligible volume of configu- 
ration space, they will not be sampled at all. For this 
reason, it is imperative to check this method for small 
systems for which all distinct basins can be identified. 

To test the method, we consider N disks in a square 
box of length L with periodic boundary conditions. 
Disks i and j interact via a "harmonic" repulsion Vij = 
e (1 — /aij)'^ /2 when the distance between their cen- 
ters of mass, rij is smaller than the sum of their radii, 
(Tjj = ((Ti -|- (7j) /2, and zero otherwise. In order to avoid 
crystallization, we use a 50 : 50 binary mixture of disks. 
The diameter ratio of the large disks to the small ones 
is 1.4. We choose units where the length of the simula- 
tion box is L = 1 and the characteristic energy of the 




b tot 



FIG. 1: (color online) Probability of finding a given minimum 
calculated in two ways: from direct enumeration, Pd, and 
from MC calculations of the basin volume relative to the total 
volume of configuration space, Vt/Vtot- Included are systems 
at packing fraction = 0.9 of A'' = 8 (black circles), 10 (red 
squares), 12 (blue diamonds), 14 (green upward triangles), 
and 16 (orange pluses) particles. For A*' = 8, 10, and 12, 
all distinct states are shown, while for A'' = 14 and 16 only 
the first 1000 states are shown. The dashed black line is 
Pd = Vb/Vtot- Inset: the volumes of all distinct basins for 
= 8, calculated by steepest descent (vf ^, red triangles) and 
conjugate gradient {v^'^ , black squares) compared to volumes 
Vb calculated by the L-BFGS algorithm. The dashed black 
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interaction is e = 1. For this system, the total volume of 
configuration space occupied by distinct basins is 



-dN 



[{N/2)\ 



(5) 



where [(A^/2)!] accounts for disk indistinguishability. 

The direct calculation of the integral on the right hand 
side of Eq. ([3]) is computationally expensive because the 
acceptance step of every MC move requires an energy 
minimization (to see if the system has left the original 
basin). Otherwise, the calculations are exactly as in 
Ref. [l^. In what follows, we use Gauss-Lobatto quadra- 
ture to evaluate Eq. (O, changing variables so that the 
integrand varies only weakly over the integration inter- 
val to improve accuracy (see [31). We verified that the 
integrand in the Gauss-Lobatto quadrature indeed varies 
smoothly with increasing force constant k of the har- 
monic spring. 

As the potential energy has to be minimized at ev- 
ery step, the efficiency of the energy minimization rou- 
tine becomes important. From any given starting point, 
the routiire should find the minimum corresponding to 
a steepest-descent (SD) search. Only the SD algorithm 
itself is guaranteed to do that, but this algorithm is not 
efficient at finding the minimum. In what follows, we 
make use of the L-BFGS minimization routine |18| as 
it is much (an order of magnitude) faster. We find that 
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FIG. 3: (color online) Configurational entropy Sc/ks as a 
function (a) system size N a.t tf) = 0.9, and (b) packing fraction 
(fi at N = 16. The red line in (a) is the linear fit to the data: 
Sc = 0.83N - 2.48. 



FIG. 2: (color online) Cumulative distribution of the basin 
volume, ^i-^) for binary mixtures with A'^ = 10 (black cir- 
cles), 12 (red squares), 14 (blue diamonds), and 16 (purple 
triangles), all at (p — 0.9. The orange solid curve shows the 
quasi log-normal fit to the = 16 data according to Eq. (|6]) 
with a = 0.23, b = 0.60, and c = 1.04. 



the L-BFGS, conjugate gradient (CG) and SD algorithms 
yield very similar results for basin volumes and volume 
distributions for iV = 8, as shown in the inset to Fig. [T] 
However, in general it may be safer to use SD, in spite of 
its higher computational cost. 

The first step in the computation of a basin volume is 
to find a potential energy minimum. To do this, we gen- 
erate a random point in the configuration space of the 
system under study (a diV-dimensional hypercube for a 
system of N spheres in d spatial dimensions). Starting 
from this initial coordinate, the potential energy of the 
system is minimized to find the coordinate Rq that cor- 
responds to the (local) potential minimum Since the 
probability of sampling a given minimum is proportional 
to the volume of its "catchment basin" , we can deduce the 
volumes of the individual basins from the frequency with 
which they are sampled, for systems sufficiently small 
so that all basins can be sampled in a simulation. Thus, 
this brute-force approach can be used to validate the free- 
energy based volume calculation. 

We used the two approaches mentioned above to com- 
pute the number of distinct catchment basins, flc, of the 
binary disk mixture at a packing fraction (j) = 0.9 and 
system sizes N G [8,16]. For these small systems, we 
can find effectively all distinct states by sampling up to 
Nt = 10* uncorrelated initial configurations 15|. During 
the runs, we keep track of ns{nt), the number of distinct 
basins sampled after nt randomly chosen initial configu- 
rations. As shown in Ref. [isj . Us saturates for large rit, 
suggesting that we have found all distinct basins or, more 
precisely: the combined volume of all basins not sampled 



is less than 0(n^^). The fractional volume occupied by 
an individual basin i is then given by Pd{i) = n{i)/Nt, 
where n{i) denotes the number of times that we have 
sampled the same basin i after Nt trials. 



For each distinct basin, we also calculate the basin vol- 
ume Vb{i) using the free-energy method described above. 
The fractional volume occupied by distinct basin i is 
given by Vb{i)/Vtot, where Vtot is given by Eq. ©. 

Fig. [T] shows the correlation between Pd and Vb/Vtot 
for each of the distinct basins i obtained from the direct 
enumeration. The dashed line is not a fit but corresponds 
to the relation = Vb/Vtot- Thus, Fig. [T] shows that 
the free-energy calculation of the basin volumes works 
very well, even though the shapes of the high-dimensional 
basins are very complicated. 

It is straightforward to calculate the average basin vol- 
ume (vb) when all the distinct basins are known. But for 
larger systems for which only a small subset of basins can 
be identified, {vb) can be calculated only if the distribu- 
tion of basin volumes, P{vb), scales in a known fashion 
with system size. Fig. [2] shows the cumulative distribu- 
tion of the basin volumes For larger N, the 

cumulative distribution P{vb) is well represented by 



{erf [ahi{x)+b] + ly 



(6) 



where x = '^^^ while a, b, and c are adjustable param- 
eters. As N increases, a and c decrease slightly, while b 
increases. Specifically, c ^ 1, suggesting that the distri- 
bution P(Df,)bccomes log-normal for larger systems. 

This result is perhaps not surprising if one expects the 
distribution of the entropy of states within a basin (the 
logarithm of the basin volume) to be Gaussian in the 
thermodynamic limit. If this is indeed the case, then one 
can compute the average basin volume (and hence the 
total number of distinct basins) from a simulation that 
samples only a fraction of all basins. 

Once I{vb) has been obtained for a given system size, 
the configurational entropy Sc follows, using Eq. Q. 
Fig. [21^ a) shows that the configurational entropy Sc = 
/cs In is extensive, i.e. it scales linearly with N. This 
is expected for large systems d, [l^ but not necessarily 
for the sizes studied here. Fig. ^h) shows the variation 
of Sc with the packing fraction cf). The number of distinct 
states increases as cj) decreases, as expected. Note that 
Sc{4>) is the configurational entropy of distinct jammed 
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FIG. 4: Entropy of distinct minima at packing fraction (j) 
and pressure p, Sc{<i>,p) /ks = lnil{<f),p), where fl{cf),p) is the 
density of states between p and p + dp. Here, Sc{4'-iP) is shown 
for systems of = 16 particles at </!> = 0.82 (solid), 0.83 
(dashed), 0.84 (dot-dashed) and 0.85 (dotted). Sc reaches a 
well-defined value as p — ^ 0. 

energy minima in soft sphere packings, or equivalently, 
the entropy of distinct mechanically-stable packings. 

However, the entropy of distinct mechanically-stable 
packings of hard spheres, So{4>)i is not the same as that 
for soft spheres iS'c(0). To obtain the former quantity, 
we must look only at soft-sphere packings that are at 
the jamming threshold {i.e. at zero pressure, p — 0) 
at each packing fraction Fortunately, we can calcu- 
late this directly from the sampled soft-sphere minima 
without introducing a protocol for bringing the system 
to p = that might bias the weightings of the result- 
ing states [l3|. We calculate the distribution P{(j>,p) 
of basins whose minima have pressure p sd (j) and use 



the average basin volume, (vb), to obtain the density 
of states of distinct energy minima, fl{(l),p), with pres- 
sures between p and p + dp, via Eq. (j4l). The en- 
tropy of distinct jammed hard-sphere packings is then 
So(,(b) = SM,p = 0) = kB lnr!(0,p = 0). 

Fig. 13] shows that S'o((/)) increases with decreasing 
over the range studied. However, we also know that Sq 
must vanish at sufficiently small </>. Thus, Sq must have 
a maximum, in agreement with earlier estimates [Tsj and 
theoretical predictions Q. It would be interesting to ex- 
plore the connection between this maximum and the ran- 
dom close- packing density in large systems [l^] . 

In summary, the free-energy method proposed here al- 
lows us to compute the volume of individual basins in 
the energy landscape of a many-particle system. This, 
in itself, is an extremely useful result. We also find that 
from the distribution of basin volumes we can obtain the 
number of distinct energy minima (the number of distinct 
jammed packings). Here, we have tested our method for 
small systems where all basins can be identified by brute 
force, but our method can be applied to far larger sys- 
tems, where direct enumeration is impossible. In prac- 
tice, the reliability of this approach depends strongly on 
the existence of a universal form for the functional form of 
the distribution basin volumes. Further tests are needed, 
but our results suggest that a log-normal form may be 
appropriate for larger system sizes. 
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